﻿#pragma once

#include <xy/math/matrix.h>

namespace xy
{

namespace solver
{

class linear
{
  public:
    enum class Status
    {
        Uninitialized,
        Success,
        Failure,
        Singular,
    };

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear(xy::math::matxd &&A);

    /**
     * @brief 获得状态
     *
     * @return Status
     */
    Status status() const;

    /**
     * @brief 获得状态
     *
     * @return Status&
     */
    Status &status();

    /**
     * @brief 移动矩阵
     *
     * @return xy::math::matxd&&
     */
    xy::math::matxd &&move_matrix();

    /**
     * @brief 计算接口
     *
     */
    virtual void perform();

    /**
     * @brief 求解接口
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    virtual xy::math::matxd solve(xy::math::matxd &&b) = 0;

  protected:
    xy::math::matxd A;
    Status m_status;
};

class linear_l : public linear
{
  public:
    linear_l(linear_l &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_l(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_l(xy::math::matxd &&A);

    /**
     * @brief 求解线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t row, std::size_t col);
};

class linear_u : public linear
{
  public:
    linear_u(linear_u &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_u(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_u(xy::math::matxd &&A);

    /**
     * @brief 求解线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t row, std::size_t col);
};

class linear_lu : public linear
{
  public:
    linear_lu(linear_lu &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_lu(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_lu(xy::math::matxd &&A);

    /**
     * @brief LU 分解，返回状态
     *
     * @return void
     */
    void perform() override;

    /**
     * @brief 求解线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t row, std::size_t col);
};

class linear_ptlu : public linear
{
  public:
    linear_ptlu(linear_ptlu &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_ptlu(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_ptlu(xy::math::matxd &&A);

    /**
     * @brief LU 分解，返回状态
     *
     * @return void
     */
    void perform() override;

    /**
     * @brief 求解线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t row, std::size_t col);

  protected:
    std::vector<std::size_t> m_perm;
};

class linear_llt : public linear
{
  public:
    linear_llt(linear_llt &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_llt(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_llt(xy::math::matxd &&A);

    /**
     * @brief LLT 分解，返回状态
     *
     * @return void
     */
    void perform() override;

    /**
     * @brief 求解线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t row, std::size_t col);
};

class linear_ldlt : public linear
{
  public:
    linear_ldlt(linear_ldlt &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_ldlt(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_ldlt(xy::math::matxd &&A);

    /**
     * @brief LLT 分解，返回状态
     *
     * @return void
     */
    void perform() override;

    /**
     * @brief 求解线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t row, std::size_t col);
};

class linear_thomas : public linear
{
  public:
    linear_thomas(linear_thomas &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_thomas(const xy::math::diagmat3d &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_thomas(xy::math::diagmat3d &&A);

    /**
     * @brief 三对角分解，返回状态
     *
     * @return void
     */
    void perform() override;

    /**
     * @brief 求解线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::tuple<xy::math::diagmat3d, xy::math::matxd, xy::math::matxd>
     */
    static std::tuple<xy::math::diagmat3d, xy::math::matxd, xy::math::matxd> example(std::size_t row, std::size_t col);
};

class linear_blind : public linear
{
  public:
    linear_blind(linear_blind &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_blind(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_blind(xy::math::matxd &&A);

    /**
     * @brief 盲人爬山法分解，返回状态
     *
     * @return void
     */
    void perform() override;

    /**
     * @brief 计算矩阵范数 ||A^-1||_infty = ||b||_1
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, A^-1, x0
     *
     * @param[in] n
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t n);
};

class linear_optimal : public linear
{
  public:
    linear_optimal(linear_optimal &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     * @param[in] x
     */
    linear_optimal(const xy::math::matxd &A, const xy::math::matxd &x);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     * @param[in] x
     */
    linear_optimal(xy::math::matxd &&A, xy::math::matxd &&x);

    /**
     * @brief 优化线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t row, std::size_t col);

  protected:
    xy::math::matxd x;
};

class linear_jacobi : public linear
{
  public:
    linear_jacobi(linear_jacobi &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_jacobi(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_jacobi(xy::math::matxd &&A);

    /**
     * @brief 求解线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t row, std::size_t col);

  public:
    double epsilon = xy::math::numeric::epsilon;
};

class linear_gauss_seidel : public linear
{
  public:
    linear_gauss_seidel(linear_gauss_seidel &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_gauss_seidel(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_gauss_seidel(xy::math::matxd &&A);

    /**
     * @brief 求解线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t row, std::size_t col);

  public:
    double epsilon = xy::math::numeric::epsilon;
};

class linear_sor : public linear
{
  public:
    linear_sor(linear_sor &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_sor(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_sor(xy::math::matxd &&A);

    /**
     * @brief 求解线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t row, std::size_t col);

  public:
    double omega = 1.2;
    double epsilon = xy::math::numeric::epsilon;
};

class linear_qr : public linear
{
  public:
    linear_qr(linear_qr &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_qr(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_qr(xy::math::matxd &&A);

    /**
     * @brief QR 分解，返回状态
     *
     * @return void
     */
    void perform() override;

    /**
     * @brief 求解线性方程组的解，b 以行为输入，返回 x 以行为输出
     *
     * @param[in] b
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&b) override;

    /**
     * @brief 示例接口 A, x, b : A -> col x col * random, x, b -> row x col
     *
     * @param[in] row
     * @param[in] col
     * @return std::array<xy::math::matxd, 3>
     */
    static std::array<xy::math::matxd, 3> example(std::size_t row, std::size_t col);

  protected:
    xy::math::vecxd beta;
};

class linear_eigen : public linear
{
  public:
    linear_eigen(linear_eigen &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_eigen(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_eigen(xy::math::matxd &&A);

    /**
     * @brief 求解对应特征值的特征向量并返回
     *
     * @param[in] eigenvalue
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&eigenvalue) override;
};

class linear_sym_eigen : public linear
{
  public:
    linear_sym_eigen(linear_sym_eigen &&) = delete;

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_sym_eigen(const xy::math::matxd &A);

    /**
     * @brief 构造函数
     *
     * @param[in] A
     */
    linear_sym_eigen(xy::math::matxd &&A);

    /**
     * @brief 求解对称矩阵的特征值和特征向量，返回特征向量构成的正交阵
     *
     * @param[in] sigma 关值
     * @return xy::math::matxd
     */
    xy::math::matxd solve(xy::math::matxd &&sigma) override;

  protected:
    /**
     * @brief 沿着 p,q 分量旋转得到 JT*A*J, Q 累计 J
     *
     * @param[in] Q
     * @param[in] p
     * @param[in] q
     */
    xy::math::matxd rotate(xy::math::matxd &&Q, std::size_t p, std::size_t q);

  protected:
    double sigma;
};

} // namespace solver

} // namespace xy